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Abstract 

Cells of Voronoi diagrams in two dimensions are usually considered as having edges 
of zero width. However, this is not the case in several experimental situations in 
which the thickness of the edges of the cells is relatively large. In this paper, the 
concept of a thick Voronoi tessellation, that is with edges of non-zero width, is 
introduced and the statistics of cell areas, as thickness changes, are analyzed. 
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1 Introduction 



Voronoi tessellations provide a powerful method for subdividing space in ran- 
dom partitions and have been used in a such diverse fields as statistical me- 
chanics 0, quantum field theory astrophysics jij, structure of matter 
0, @, 0, biology Q and medicine @. 



Algorithms to generate Voronoi diagrams can be found, among others, in [10 



111 ]. [121 ] [131 ]; usually Voronoi partitions are obtained by considering a uniform 
spatial distribution of centers, so that the probability that there are n points 
in a given domain of space obeys a Poisson distribution, of constant intensity 
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A 12J. Non Poissonian distribution of seeds have been used in 14|, in an 



astrophysical framework, and examples of area distribution for this case can 
be found in 1151. 



Usually, when generating Voronoi diagrams, the thickness of the edges is not 
taken into account, and, hence, its influence on the size of the cells is not con- 
sidered. However, in nature, configurations occur that approximate Voronoi 
tessellations in which the width of edges is not negligible; examples can be 
found in the generation of diffraction patterns 16] , in the compaction of gran- 
ular matter 17 [and crystallization of granular fluids 18], and in animal coat 



formation [19|, [20 



In the next section probability distribution used to fit areas of 2D Poissonian 
Voronoi tessellation will be briefly reviewed; next the statistics of Poissonian 
thick Voronoi diagrams, that is with an edge of non-zero width, will be studied 
in the two-dimensional case and, in particular, the dependence of the mean 
and variance of the area distribution on the thickness of edges will be analyzed. 
Furthermore, a probability density function will be considered, adapted from 



that proposed in [2l|], to fit histograms of cell areas. 



2 Probability distributions 



Consider a. 2D Poisson Voronoi Diagram (PVD for short) and let ao be the 
area of its cells. 



No exact solution is known for the distribution of 2D Voronoi cells, however 
there exist several PDF, based of the Gamma distribution, that provide an 
approximate solution, with different degrees of "goodness": the most general 
way is via a 3 parameter generalized Gamma distribution 



G(x; a, b, c) 



r(f) 



221 ]. 13]; here x = ao/(ao), where (ao) is the area average. The previous PDF 
can be simplified inserting a = 1 and the following 2 parameters Gamma PDF 
is obtained 



g(x;b,c) 



b c x c x e bx 



r(c) 



(2) 



[lil,l23] 
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Note that the most recent numerical estimates of the parameters c and 5, in 
the 2D case, show that they are very close, namely 6=3.52418 and c= 3.52440 
231 ] . This suggests the use of a one-parameter Gamma distribution, 



Tic) 



-x 



c-l 



exp( 



-ex) 



(3) 



such a distribution has been used by Kiang in his seminal work on Voronoi 
diagrams 24j. Note, however, that in 24j it is has been assumed c = 2d where 
d is just the dimensionality of the cells, so in the present case case c = 4. 



Finally in [2l|, a simpler distribution has been proposed with no free param- 
eters, that can be obtained from (3) by setting c = (3d + l)/2, that is c = 3.5 
in 2D: 



f(x; d) = Cx~ exp (-(3d + l)x) 



(4) 



where 



C 



3d+l 



(3<f+l)/2 



2 J 



(5) 



A detailed comparison between f(x; d) and G(x; a, b, c) can be found in [21 
The usefulness of these distributions in physical problems is determined by 
the trade off between the accuracy with which they fit the data and their 
complexity. We have computed 5 = \f(x; d) — g(x; b, c)| the absolute value of 
the difference between f(x; d) and g(x; b, c), for the 2D and the resulting plot 
is shown in Fig. 1: the maximum value of 5 is 0.0037. Obviously, the PDF g 
can give a better fit of the data, however / is simpler and we think that, at 
least for the purposes of this note, provides a good enough approximation. 

For future reference we give here the variance and mode of (4) which are, 
respectively, 

° 2 - ii+i < 6 > 



and 



3d 



3d+l 



(7) 
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Fig. 1. The absolute value of the difference between f(x;d) and g(x;b,c) as a 
function of x when d = 2, 5=3.52418 and c= 3.52440. 

3 Two-dimensional thick Voronoi diagrams 



Let s be the thickness of edges of 2D cells in a Voronoi tessellation, denote 
with a(s) the cells area and with ao the area when s = 0. 

The analysis of area size as a function of s can be made independent of the 
area A of the 2D domain A in which the Voronoi polygons are generated, by 
introducing a dimensionless parameter 

s s 



where n is the number of seeds of the diagram. 

An example of thick Voronoi diagrams, with p = 0.2, is presented in Fig. 
2 where, for illustrative purposes, just n = 150 centers have been used; all 
simulations presented in the following have been carried out with n = 4 • 10 4 
centers. 



Simulations were run on a LINUX -2.66GHz processor: Poisson Voronoi tes- 
sellation were generated by sampling independently the coordinates along X 
and Y axes from a uniform distribution by means of the subroutine RAN2 



described in [25(. In order to minimize boundary effects introduced by cells 
crossing the boundary of the domain A, a square B is defined larger by a 
factor 1.5 and containing A. The seeds are placed in the whole B domain; 
furthermore only cells that do not cross the boundary of A are considered, see 
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Fig. 2. Thick Voronoi diagrams in 2D with p = 0.2. The selected region comprises 
150 random seeds marked by a point. Here, for simplicity (ao) = 1. 

Figure 2. 



Further information on the code used here can be found in |26|]. The CPU 
time running time was 1.42 s for a seed. 

It is obvious that increasing values of p make more likely the occurrence of 
cells completely covered by edges and indeed, from Fig. 2, it can be seen that 
some of the smallest cells are completely covered. Moreover it is also apparent 
for p = 0.2 a relatively large area of the cell is occupied by edges and this is 
confirmed by the results presented in the sequel. We have then restricted our 
analysis to < p < 0.2. 

In general the area a can be considered to be obtained from ao = a(0) by the 
action of a mapping K, so that a(p) = K(a ,p) rescales the cell size. Then K 
can be seen as a nonlinear scaling operator and can be given in a form akin 
to that used for the group of scaling [27l |. namely 



a(p) = k(a ,p)a . 



(9) 



Note that A; is a decreasing function of ao, in that large cells are relatively less 
affected than small ones by the occurrence of an edge of width s. 



Since cells are irregular polygons it is difficult to compute an explicit form 
of a and k since it depends on the shape of the cell, however some general 
properties are readily apparent, which suffice for our purposes. The area a 
can not contain powers of s, and hence of p, larger than 2, for reasons of 
dimensional consistency, and that holds for k too (compare Eqs. (9)); it is also 
obvious that 

fc(a ,0) = l, dk(a o ,p)/dp\ p=o <0. (10) 

Then k can be written as 

k(a ,p) = 1 - h(a )p + ^g(a )p 2 , (11) 

where 

h(a ) =| dk(a ,p)/dp \ p=0 \, g(a ) = d 2 k(a , p)/dp 2 \ p=0 , 

and the area a(p) is 

«(p) =ao(l- h(a )p + i#(a )p 2 ) , (12) 

with the understanding that a(p) is set to if Eq.(12) yields a value less then 
0. 

Since p is small, only the linear term in (11) needs to be considered, then 
a(p) = a (l - h(a )p + ^g(a )p 2 ^j « a - a h(a )p, (13) 

from which 

(a(p)) w (ao) - (a h(a ))p. (14) 

A linear decrease of the average area is also the outcome of simulations, as 
shown in Fig. 3 where a comparison with a linear fit is presented: values of 
(a ) and (a h(a )) have been computed by standard fitting procedures (least 
squares method) and are reported in Table 1. 

It should be observed that the area decreases quite sharply even for values of 
p which allow the linear approximation of Eq. (13); in particular for p = 0.2, 
the average area is (a(p)) ~ 0.5(a ). The trend is linear up to p = 0.25 
approximately, for this value (a(p)) ~ 0.3(ao). For illustrative purposes Fig. 4 
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Fig. 3. Area average versus thickness. The stars represent the results of the simu- 
lations and the full line is given by Eq. (14). 

Table 1 

Coefficients of the linear fit. 



(a ) (h(a )a } a 2 C 



0.19 0.55 0.01 0.016 



shows the trend of (a(p)) in the interval < p < 0.5, fitted with the equation 
obtained by averaging the terms of Eq. (12), that is 

(a(p)) = (a ) - (a h{a ))p + ^(a g{a ))p 2 . (15) 
The variance a 2 can then be calculated: 



o"a = (( a - ( a )) 2 ) = ([K - a h(a )p) - ((a - a h(a ))} p} 2 ) , (16) 



from which it is straightforward to obtain: 



°"a = ((°o - (a )) 2 ) - 2 ((h(a )a - (h(a )a })(a - (a ))) p 




(17) 



and, by making use again of the linear approximation, 

o 2 a w <J(a - (a )) 2 ) - 2p ((h(a )a - (/i(a )a ))(a - (a ))) . 



(18) 




Fig. 4. Area average versus thickness. The stars represent the results of the 
simulations and the full line Eq. (15) with (a = 0.2) , (/i(ao)ao) = 0.68 and 
l/2<a s(ao)> = 0.67. 



Set 



{(h(a )a - (h(a )a ))(a - (ao))) = C(h(a )a , a ), 



then Eq.(18) can be written as 



a a ~ a a ~ 2C(h(a )a , a )p. 



(19) 



Now C(h(ao)aQ, ao) > 0, because it is the covariance of h(ao)ao with a and, 
by definition, both a and /i(a ) are positive; therefore, if the approximation 
of Eq. (18) holds, a 2 a must fall off linearly. Comparison between simulations 
and the results of Eq. (19) are shown in Fig. 5 and it is clear from the figure 
that the effect of a quadratic term on the fit is negligible; indeed we have 
verified that ((/i(ao)ao — (M a o)°o)) 2 ) is about an order of magnitude smaller 
than C(h(ao)ao,ao). Numerical values of a^ Q and C(h(ao)ao,ao) are shown in 
Table 1. 



4 Fitting area distributions of thick Voronoi diagrams 



We consider now the probability density function (PDF) of a; as noted earlier 
PDFs of Voronoi cell size are commonly expressed in terms of a standardized 
variable x, which, in the present case, takes the form x = a/ < a >, so that, 
obviously, (x) = 1. For future reference, we compute the variance <j 2 x of x, 
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Fig. 5. Variance of cells area versus thickness. The stars represent the results of the 
simulations and the line the theoretical variance, as given by formula (19). 



which is given by 



a 



(20) 



from which, making use of Eqs. (14) and (19), 



- 2C(h(a )a , a )p 



(a - a h(a )pY 



(21) 



Values of a 2 x provided by Eq. (21) are in good agreement with the results of 
simulations as shown in Fig. 6. 

To investigate the distribution of cell areas, we have used the probability 
density function (4) : however, since a 2 x increases with p, from Eq. (6) it is 
clear that, in order to use the PDF / to fit histograms of simulated data, d 
must be considered to be a variable parameter which decreases for increasing 
p. From Eqs. (6) and (21), it is straightforward to derive a formula for d: 



2 (q - a h(ao)p) 2 
a l ~ 2C(h(a )a ,a )p 



(22) 



Empirical values of d have been found by the method of matching moments 
and are shown in Fig. 7, together with the fit provided by Eq. (22). We have 
generated histograms of area distribution for different values of p and have 
used the PDF / with the corresponding parameter d derived with Eq. (22): 
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Fig. 6. Variance of the standardized variable x versus thickness. The stars represent 
the results of the simulations and the full line reports the theoretical variance as 
given by formula (21). 
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Fig. 7. The dimension d that models the Voronoi cell standardized area-distribution 
in 2D as function of the adimensional parameter p (dotted line), see formula (22) 
and simulated points (stars) . 

statistical tests show that xl increases with p and that xl < 1-53 up to 
p = 0.04 (see Fig. 8), thus the fit is adequate only for very small values of p. 

Even though the PDF /, with d given by Eq. (22), gives good results only for 
small p values, it can be used to predict, at least qualitatively, the change in 
shape of the empirical distribution. The decrease of the parameter d implies 
a shift of the mode a m toward zero (see Eq.( 7)) and that occurs also in the 
histograms generated by the simulations. 
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Fig. 8. Histogram (step-diagram) of the Voronoi normalized thick area distribution 
in 2D with a superposition of the gamma PDF derived by Eq. (4), with d given 
by Eq. (22). The number of seeds and bins are 40000 and 100 respectively: here 
d = l.n ,p = 0.04, X 2 = 151 ^ = 1-53. 




The result is the appearance of a PDF decreasing monotonically with a, as an 
example see Figure 9, where is clear that lack of agreement with the modified 
PDF /. 
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5 Conclusions 



Voronoi Diagrams in 2D are usually generated as irregular polygons whose 
edges, in principle, have zero thickness: however in several experimental situ- 
ations, Voronoi cells appear to have edges of relatively large width. 

Clearly the emergence of thick edges is related to the formation of configura- 
tions representing approximate Voronoi diagrams, as results of chemical and 
phy sical mechanisms, are in general quite complex. 

Consider, for instance, pattern formation in certain animals coats (e.g. giraffe) 
by reaction diffusion processes. Here two-dimensional Voronoi diagrams are 
generated by an activator a, diffusing from randomly placed point sources, 
which switches on the production of melanin, this switch being controlled 



by a threshold 9 19j; in a more complex model [20| melanin production is 



modulated by the concentration of a substrate s. Thickness of edges is then 



determined by the value of 9 [19(, or by the abundance and removal rate of 



the substrate 20 . 



Voronoi cells can also be generated when an homogeneous medium is occupied 
by domains emerging from random placed seeds with the same isotropic growth 
rate 28|,|29|: such is the case of crystals 28|,|29[ or bubbles in volcanic eruptions 



301 ] . In this case the width of the edges can be determined by the relation 
between the amount of growth and the size of the domain where it takes 
place. 

Edges width should be constant, at least approximately, when processes lead- 
ing to the formation of Voronoi cells are symmetric, whereas if symmetry 
breaks down edges of different thickness must be expected. Consider again 
animal coat formation: if the threshold 9 is not constant over the domain 
where cells are formed edges of different width will emerge. The same effect 
can result if the constant value of 9 is replaced by a probability distribution 
p(9). 

Here a general method has been presented to compute the statistics of cell 
areas as the thickness parameter p varies: here, for each value of p, edges have 
the same width. Theoretical computations as well as results of simulations 
show that the mean area (a) and variance a\ fall off linearly for p increasing 
in the interval [0,0.2]: in particular the mean area shows a marked decrease 
and at p = 0.2 is reduced to 50% of its original value. 

We have tried to fit the simulated distribution of the standardized variable 



x = a/ (a) with the PDF presented in [2l|, by using d as a free parameter, but 



such a fit holds only for p < 0.07; for larger values of p simulations show that 
the mode shifts close to zero more rapidly than predicted by equation Eq. (7). 
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Finally, it should be noted that the formation of thick edges can be seen as 
a particular example of a process by which cells are eroded. The approach 
presented here, however, is general enough to be readily adapted to analyze 
the statistics of cell areas for different cases of cell erosion, for instance as 
in diffusion-limited aggregation of Voronoi diagrams jsj , in that every erosion 
operator k must have the form given by Eq. (11) and the terms h and g can 
be determined from the data, experimental or simulated. 
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